grid_5km_v1 <- read_sf(here('spatial_data', 'master_5km_grid_tmer.shp'))
grid_5km_v2 <- read_sf(here('spatial_data', 'fivekm_grid_polys_shore_lamb.shp'))14_blh_new-5km-shapefile
I ran the monthly maps code and noticed a couple differences, all are reasonable:
Grid cells are trimmed to the coastline. Great, much prettier maps.
New BC subdirectory. Makes sense for a couple grid cells at the northern end of Washington, and will go away once I filter out the Puget Sound.
A couple of grid cells moved from northern end of OR map to southern end of WA map. Seems fine, not a huge impact on the output, and the data doesn’t seem to be lost, just moved along the border.
Not related to the change in shape files, but related to the map visualizations - We still need to decide if it makes more sense to facet based on fishing effort (where the VMS ping is, how it’s done now) vs. landings (where the ticket is landed, how it could be done). I think the answer should probably be similar to our decision to visualize by whale season, what makes the most sense from a management perspective? What regulations are the fishers held to, is it based on area fished or area landed?
I want to get a sense of the unique values and combinations of values for the new features, compared to the old features.
I’m going to load the two versions and compare.
What dimensions?
dim(grid_5km_v1)[1] 76260 10dim(grid_5km_v2)[1] 72721 41How many unique grid cells are there?
n_distinct(grid_5km_v1$GRID5KM_ID)[1] 76260n_distinct(grid_5km_v2$GRID5KM_ID)[1] 72721There are more grid cells in the older version, which makes me wonder if the extent is slightly broader. There are 31 more columns in the newer version, which makes sense because Blake made this to include more features I asked for.
What features do they each have?
colnames(grid_5km_v1) [1] "GRID5KM_ID" "centro_lat" "centro_lon" "STATE" "WM_NGDC_M"
[6] "EEZ" "WATER" "Centroid_y" "Centroid_x" "geometry" colnames(grid_5km_v2) [1] "GRID5KM_ID" "ORIG_AREA" "AREA" "Shoreline" "WM_NGDC_m"
[6] "NGDC_count" "NGDC_med_m" "NGDC_SD" "NGDC_min_m" "NGDC_max_m"
[11] "NGDC_var" "STATE" "EEZ" "Water_zone" "Salish_Sea"
[16] "RAMP_area" "RAMP_zone" "WDFW_WSMA" "EEZ_MRGID" "Derville"
[21] "centro_lat" "centro_lon" "Centroid_y" "Centroid_x" "LL_EAST"
[26] "UL_EAST" "UR_EAST" "LR_EAST" "LL_NORTH" "UL_NORTH"
[31] "UR_NORTH" "LR_NORTH" "LL_LATITUD" "UL_LATITUD" "UR_LATITUD"
[36] "LR_LATITUD" "LL_LONGITU" "UL_LONGITU" "UR_LONGITU" "LR_LONGITU"
[41] "geometry" # what is the same?
intersect(colnames(grid_5km_v1), colnames(grid_5km_v2))[1] "GRID5KM_ID" "centro_lat" "centro_lon" "STATE" "EEZ"
[6] "Centroid_y" "Centroid_x" "geometry" # what is unique to V1?
setdiff(colnames(grid_5km_v1), colnames(grid_5km_v2))[1] "WM_NGDC_M" "WATER" # what is unique to V2?
setdiff(colnames(grid_5km_v2), colnames(grid_5km_v1)) [1] "ORIG_AREA" "AREA" "Shoreline" "WM_NGDC_m" "NGDC_count"
[6] "NGDC_med_m" "NGDC_SD" "NGDC_min_m" "NGDC_max_m" "NGDC_var"
[11] "Water_zone" "Salish_Sea" "RAMP_area" "RAMP_zone" "WDFW_WSMA"
[16] "EEZ_MRGID" "Derville" "LL_EAST" "UL_EAST" "UR_EAST"
[21] "LR_EAST" "LL_NORTH" "UL_NORTH" "UR_NORTH" "LR_NORTH"
[26] "LL_LATITUD" "UL_LATITUD" "UR_LATITUD" "LR_LATITUD" "LL_LONGITU"
[31] "UL_LONGITU" "UR_LONGITU" "LR_LONGITU"The newer version has more information about:
Area. The new version includes the area assuming the grid cell is 5km x 5km, and the actual area trimmed to exclude the coastline.
Shoreline. Binary for whether a grid cell intersects the shoreline.
Depth. The older version only had the weighted mean depth of the grid cell, the new version also includes the median, standard deviation, variance, minimum and maximum.
Water zone. Does the grid cell fall within <3, 3-200, >200 nautical miles from shoreline.
Salish sea. Does the grid cell fall within the Puget Sound or not.
Management areas. CDFW RAMP areas and zones for California, Derville zones for Oregon, and Marine Fish and Shellfish Manage Catch Reporting areas for Washington.
Corner coordinates. Since I want to be able to join this map to Chris Free’s dataframe which has latitude over time, Blake added the upper and lower left and right corner locations of each grid cell.
This is all based on the metadata file that Blake created, which is super helpful, since that detail is not knowable from column names alone.
What geometry type?
class(grid_5km_v1$geometry)[1] "sfc_POLYGON" "sfc" class(grid_5km_v2$geometry)[1] "sfc_MULTIPOLYGON" "sfc" Version 1 is polygons. Version 2 is multipolygons, because it trims the grid cells to the coast, and the shore can sometimes split a grid cell into multiple polygons.
What coordinate system?
st_crs(grid_5km_v1)Coordinate Reference System:
User input: WGS_1984_Transverse_Mercator
wkt:
PROJCRS["WGS_1984_Transverse_Mercator",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",31.96,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-121.6,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",390000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]st_crs(grid_5km_v2)Coordinate Reference System:
User input: CA_Curr_Lamb_Azi_Equal_Area
wkt:
PROJCRS["CA_Curr_Lamb_Azi_Equal_Area",
BASEGEOGCRS["WGS 84",
DATUM["World Geodetic System 1984",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]],
ID["EPSG",6326]],
PRIMEM["Greenwich",0,
ANGLEUNIT["Degree",0.0174532925199433]]],
CONVERSION["unnamed",
METHOD["Lambert Azimuthal Equal Area",
ID["EPSG",9820]],
PARAMETER["Latitude of natural origin",30.5,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",-122.6,
ANGLEUNIT["Degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["False easting",1000000,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",0,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["(E)",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["(N)",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]What extent?
st_bbox(grid_5km_v1) xmin ymin xmax ymax
-616809.7 -204655.9 979184.5 1929641.3 st_bbox(grid_5km_v2) xmin ymin xmax ymax
99280 -32131 1656582 2062869 Well, they’re different, but that’s otherwise meaningless to me.
What’s it look like on a map?
background_map <- ne_states(country = 'United States of America')
bbox_v1 <- st_bbox(grid_5km_v1)
bbox_v2 <- st_bbox(grid_5km_v2)
ggplot() +
geom_sf(data = grid_5km_v1, aes(fill = STATE), lwd = 0) +
geom_sf(data = background_map) +
xlim(bbox_v1[1], bbox_v1[3]) +
ylim(bbox_v1[2], bbox_v1[4])
ggplot() +
geom_sf(data = grid_5km_v2, aes(fill = STATE), lwd = 0) +
geom_sf(data = background_map) +
xlim(bbox_v2[1], bbox_v2[3]) +
ylim(bbox_v2[2], bbox_v2[4])
What are the unique values for the new management zone variables?
unique(grid_5km_v2$RAMP_area)[1] "Southern" NA "Central" "Northern"
[5] "San Francisco"unique(grid_5km_v2$RAMP_zone)[1] NA "Zone 6" "Zone 3 & Zone 7" "Zone 5"
[5] "Zone 2" "Zone 1" "Zone 4" "Zone 4 & Zone 7"
[9] "Zone 2 & Zone 7"unique(grid_5km_v2$WDFW_WSMA) [1] NA "57" "60A-1" "61" "58A" "60A-2" "60D" "59A-2" "58B"
[10] "59A-1" "60C" "59B" "60B" unique(grid_5km_v2$Derville)[1] NA
[2] "North (North OR border to Cascade Head)"
[3] "South (Bandon to south OR border)"
[4] "Central (Cascade Head to Bandon)" What do the EEZ & management zones look like on a map?
plot_grid <- function(grid, fill_var, background) {
bbox <- st_bbox(grid)
ggplot() +
geom_sf(data = grid, aes(fill = get(fill_var)), lwd = 0) +
geom_sf(data = background) +
xlim(bbox[1], bbox[3]) +
ylim(bbox[2], bbox[4]) +
labs(fill = fill_var)
}
plot_grid(grid = grid_5km_v2, fill_var = "EEZ", background = background_map)
plot_grid(grid = grid_5km_v2, fill_var = "RAMP_area", background = background_map)
plot_grid(grid = grid_5km_v2, fill_var = "RAMP_zone", background = background_map)
plot_grid(grid = grid_5km_v2, fill_var = "Derville", background = background_map)
plot_grid(grid = grid_5km_v2, fill_var = "WDFW_WSMA", background = background_map)
How many cells intersect the shoreline?
sum(grid_5km_v2$Shoreline)[1] 1557How many cells in the Salish sea?
sum(grid_5km_v2$Salish_Sea)[1] 670Are all the grid cells from the newer version, or are any lost?
length(setdiff(grid_5km_v2$GRID5KM_ID, grid_5km_v1$GRID5KM_ID))[1] 0length(setdiff(grid_5km_v1$GRID5KM_ID, grid_5km_v2$GRID5KM_ID))[1] 3539All of them are from the older version.
Where are the old grid cells?
old_grid_cell_ids <- setdiff(grid_5km_v1$GRID5KM_ID, grid_5km_v2$GRID5KM_ID)
old_grid_5km_v2 <- grid_5km_v1 %>%
mutate(is_old = GRID5KM_ID %in% old_grid_cell_ids) %>%
filter(is_old)
old_bbox <- st_bbox(old_grid_5km_v2)
ggplot() +
geom_sf(data = old_grid_5km_v2, color = 'orange') +
geom_sf(data = background_map) +
xlim(old_bbox[1], old_bbox[3]) +
ylim(old_bbox[2], old_bbox[4])
Seems like they’re all around the edges, which I won’t be using for these maps anyways, so no problems there.
I’m going to play around with some of the map aesthetics. How to rotate California coastline?
# non-rotated version
ca_grid <- grid_5km_v2 %>% filter(STATE == "CA")
plot_grid(ca_grid, "STATE", background_map)
plot(ca_grid, max.plot = 1)
# rotated version 1
# from https://r-spatial.github.io/sf/articles/sf3.html#affine-transformations
rot <- function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)
ca_grid_rotated <- ca_grid
ca_grid_rotated$geometry <- ca_grid$geometry * rot(pi * 20/180)
plot(ca_grid_rotated, max.plot = 1)
# issue: there's no CRS when it's rotated this way, so geom_sf throws an error